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Abstract 


When the ancestors of modern Eurasians migrated out of Africa and interbred with Eurasian archaic hominins, namely, 
Neanderthals and Denisovans, DNA of archaic ancestry integrated into the genomes of anatomically modern humans. 
This process potentially accelerated adaptation to Eurasian environmental factors, including reduced ultraviolet radiation 
and increased variation in seasonal dynamics. However, whether these groups differed substantially in circadian biology 
and whether archaic introgression adaptively contributed to human chronotypes remain unknown. Here, we traced the evo- 
lution of chronotype based on genomes from archaic hominins and present-day humans. First, we inferred differences in cir- 
cadian gene sequences, splicing, and regulation between archaic hominins and modern humans. We identified 28 circadian 
genes containing variants with potential to alter splicing in archaics (e.g., CLOCK, PER2, RORB, and RORC) and 16 circadian 
genes likely divergently regulated between present-day humans and archaic hominins, including RORA. These differences 
suggest the potential for introgression to modify circadian gene expression. Testing this hypothesis, we found that intro- 
gressed variants are enriched among expression quantitative trait loci for circadian genes. Supporting the functional rele- 
vance of these regulatory effects, we found that many introgressed alleles have associations with chronotype. Strikingly, 
the strongest introgressed effects on chronotype increase morningness, consistent with adaptations to high latitude in other 
species. Finally, we identified several circadian loci with evidence of adaptive introgression or latitudinal clines in allele fre- 
quency. These findings identify differences in circadian gene regulation between modern humans and archaic hominins 
and support the contribution of introgression via coordinated effects on variation in human chronotype. 
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Significance 


Interbreeding between humans and Neanderthals created the potential for adaptive introgression as humans moved 
into environments that had been populated by Neanderthals for hundreds of thousands of years. Here, we discover 
lineage-specific genetic differences in circadian genes and their regulatory elements between humans and 
Neanderthals. We show that introgressed alleles are enriched for effects on circadian gene regulation, consistently in- 
crease the propensity for morningness in Europeans, and show evidence of adaptive introgression or associations be- 
tween latitude and frequency. These results expand our understanding of how the genomes of humans and our 
closest relatives responded to environments with different light/dark cycles and demonstrate a coordinated contribution 
of admixture to human chronotype in a direction that is consistent with adaptation to higher latitudes. 


© The Author(s) 2023. Published by Oxford University Press on behalf of Society for Molecular Biology and Evolution. 
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, 
distribution, and reproduction in any medium, provided the original work is properly cited. 
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Introduction 


All anatomically modern humans (AMHs) trace their origin 
to the African continent around 300 ka (Stringer 2016; 
Hublin et al. 2017), where environmental factors shaped 
many of their biological features. Approximately 70 ka 
(Bae et al. 2017), the ancestors of modern Eurasian AMH 
began to migrate out of Africa, where they were exposed 
to diverse new environments. In Eurasia, the novel environ- 
mental factors included greater seasonal variation in tem- 
perature and photoperiod. 

Changes in the pattern and level of light exposure have 
biological and behavioral consequences in organisms. For 
example, Drosophila melanogaster that is native to 
Europe harbors a polymorphism in “timeless,” a key gene 
in the light response of the circadian system, that follows 
a latitudinal cline in allele frequency (Sandrelli et al. 2007; 
Tauber et al. 2007). The ancestral haplotype produces a 
short TIM (S-TIM) protein that is sensitive to degradation 
by light because of its strong affinity to cryptochromes 
(CRYs), photoreceptor proteins involved in the entrainment 
of the circadian clock. Insertion of a G nucleotide in the 5’ 
coding region of the gene originated approximately 10 kyr 
in Europe and created a start codon that produces a new 
long TIM (L-TIM) isoform. The L-TIM variant has a lower af- 
finity to CRY, creating a change in photosensitivity and al- 
tering the length of the period. L-TIM flies are at a higher 
frequency in southern Europe, while S-TIM flies are more 
prevalent in northern Europe. Another example is found in pa- 
cific salmon. Chinook salmon (Oncorhynchus tshawytscha) 
populations show a latitudinal cline in the frequency and 
length of repeat motifs in the gene OtsClock 7b, strongly sug- 
gesting that this locus is under selection associated with lati- 
tude and photoperiod (O'Malley and Banks 2008; O'Malley 
et al. 2010). The evolution of circadian adaptation to diverse 
environments has also been widely studied in insects, plants 
(Michael et al. 2003; Zhang et al. 2008), and fishes, but it is 
understudied in humans. Adaptive processes could have 
helped to align human biology and chronotype to new natural 
conditions. 

Previous studies in humans found a correlation between 
latitude and chronotype (morningness vs. eveningness) 
variation (Leocadio-Miguel et al. 2017; Randler and 
Rahafar 2017; Lowden et al. 2018) and a latitudinal cline 
in some circadian allele frequencies (Dorokhov et al. 
2018; Putilov et al. 2018, 2019), highlighting the contribu- 
tion of the environment to behavior and circadian biology. 
Many human health effects are linked to the misalignment 
of chronotype (Knutson and von Schantz 2018), including 
cancer, obesity (Gyarmati et al. 2016; Papantoniou et al. 
2016, 2017; Gan et al. 2018; Shi et al. 2020; Yousef 
et al. 2020), and diabetes (Gan et al. 2015; Larcher et al. 
2015, 2016). There is also evidence of a correlation be- 
tween evening chronotype and mood disorders, most 


notably seasonal affective disorder, depression, and wor- 
sening of bipolar disorder episodes (Srinivasan et al. 
2006; Kivela et al. 2018; Taylor and Hasler 2018). Thus, 
we hypothesize that the differences in geography and en- 
vironment encountered by early AMH populations moving 
into higher latitudes created potential for circadian mis- 
alignment and health risk. 

Although AMHs arrived in Eurasia ~70 ka, other homi- 
nins (e.g., Neanderthals and Denisovans) lived there for 
more than 400 ka (Arnold et al. 2014; Meyer et al. 2014, 
2016). These archaic hominins diverged from AMHs around 
700 ka (Meyer et al. 2012; Prüfer et al. 2014, 2017; Nielsen 
et al. 2017; GOmez-Robles 2019; Mafessoni et al. 2020), 
and as a result, the ancestors of AMHs and archaic homi- 
nins evolved under different environmental conditions. 
While there was substantial variation in the latitudinal 
ranges of each group, the Eurasian hominins largely lived 
at consistently higher latitudes and, thus, were exposed 
to higher amplitude seasonal variation in photoperiods. 
Given the influence of environmental cues on circadian 
biology, we hypothesized that these separate evolutionary 
histories produced differences in circadian traits adapted to 
the distinct environments. 

When AMH migrated into Eurasia, they interbred with 
the archaic hominins that were native to the continent, ini- 
tially with Neanderthals (Green et al. 2010; Villanea and 
Schraiber 2019) around 60 ka (Sankararaman et al. 2012; 
Skoglund and Mathieson 2018) and later with Denisovans 
(Jacobs et al. 2019). Due to this, a substantial fraction 
(>40%) of the archaic variation remains in present-day 
Eurasians (Vernot and Akey 2014; Skov et al. 2020), al- 
though each human individual carries only ~2% DNA of ar- 
chaic ancestry (Vernot et al. 2016; Prüfer et al. 2017). Most 
of the archaic ancestry in AMH was subject to strong 
negative selection, but some of these introgressed alleles 
remaining in AMH populations show evidence of adapta- 
tion (Racimo et al. 2015; Gower et al. 2021). For example, 
archaic alleles have been associated with differences in 
hemoglobin levels at higher altitude in Tibetans, immune 
resistance to new pathogens, levels of skin pigmentation, 
and fat composition (Huerta-Sanchez et al. 2014; Racimo 
et al. 2015; Racimo, Gokhman, et al. 2017; Dannemann 
and Kelso 2017; Racimo, Marnetto, and Huerta-Sanchez 
2017; McArthur et al. 2021). Previous work also suggests 
that introgressed alleles could have adaptively influenced 
human chronotype. First, a phenome-wide association 
study in the UK Biobank found loci near ASB7 and EXOC6 
with introgressed variants that were significantly associated 
with self-reported sleeping patterns (Dannemann and 
Kelso 2017). One of these alleles showed a significant asso- 
ciation between frequency and latitude. Second, summar- 
izing effects genome-wide, introgressed alleles are also 
moderately enriched for heritability of chronotype com- 
pared to nonintrogressed alleles (McArthur et al. 2021). 
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These results suggest a potential role for introgressed al- 
leles in adaptation to pressures stemming from migration 
to higher latitudes. 

Motivated by the potential for a role of archaic introgres- 
sion in AMH circadian variation, we explore two related 
questions: 1) Can comparative genomic analysis identify 
differences in AMH and archaic hominin circadian biology? 
2) Do introgressed archaic alleles influence human circadian 
biology? Understanding the ancient history and evolution 
of chronotypes in humans will shed light on human adap- 
tation to high latitudes and provide context for the genetic 
basis for the modern misalignment caused by the develop- 
ment of technology and night shiftwork. 


Results 


Did Archaic Hominins and Modern Humans Diverge in 
Circadian Biology? 


Following divergence ~700 ka (Nielsen et al. 2017; Gómez- 
Robles 2019), archaic hominins and AMH were geographically 
isolated, resulting in the accumulation of lineage-specific 
genetic variation and phenotypes (fig. 1). In the next several 
sections, we evaluate the genomic evidence for divergence 
in circadian biology between archaic hominin and modern 
human genomes. 


Identifying Archaic-Hominin-Specific Circadian 
Gene Variation 


With the sequencing of several genomes of archaic homi- 
nins, we now have a growing, but incomplete, catalog of 
genetic differences specific to modern and archaic lineages. 
Following recent work (Kuhlwilm and Boeckx 2019), we de- 
fined archaic-specific variants as genomic positions where 
archaic hominins (Altai Neanderthal, Vindija Neanderthal, 
and Denisovan) all have the derived allele while, in humans, 
the derived allele is absent or present at such an extremely 
low frequency that it is likely an independent occurrence. 
We defined human-specific variants as positions where all 
individuals in the 1000 Genomes Project carry the derived 
allele and all the archaics carry the ancestral allele. 

We evaluated archaic-specific variants for their ability to 
influence proteins, splicing, and regulation of 246 circadian 
genes (Materials and Methods). The circadian genes were 
identified by a combination of literature search, expert 
knowledge, and existing annotations (supplementary 
table S1 and fig. S1, Supplementary Material online; 
Materials and Methods). The core circadian clock machinery 
is composed of a dimer between the CLOCK and ARNTL 
(BMAL1) transcription factors, which binds to E-box enhan- 
cer elements and activates the expression of the Period 
(PER1/2/3) and Cryptochrome (CRY1/2) genes (fig. 10. 
PERs and CRYs form heterodimers that inhibit the positive 
drive of the CLOCK-BMAL1 dimer on E-boxes, inhibiting 


their own transcription in a negative feedback loop. 
CLOCK-BMAL1 also drives the expression of many other 
clock-controlled genes, including nuclear receptor subfamily 
1 group D members 1 and 2 (NR1D1/2), RAR-related orphan 
receptors A and B (RORA/B), and D-Box-binding PAR BZIP 
transcription factor. ROR and REV-ERB are transcriptional 
regulators of BMAL1. CK1 binds to the PER/CRY heterodi- 
mer, phosphorylating PER and regulating its degradation. 
Similarly, FBXL3 marks CRY for degradation. Beyond the 
core clock genes, we included other upstream and down- 
stream genes that are involved in the maintenance and re- 
sponse of the clock (supplementary table S1 and fig. S1, 
Supplementary Material online). 

We identified 1,136 archaic-specific variants in circadian 
genes, promoters, and distal candidate cis-regulatory ele- 
ments (cCREs). The circadian genes with the most archaic- 
specific variants are CLDN4, NAMPT, LRPPRC, ATF4, and 
AHCY (125, 112, 110, 104, and 102, respectively) 
(supplementary table S2, Supplementary Material online). 


Fixed Human- and Archaic-Specific Variants Are 
Enriched in Circadian Genes and Associated 
Regulatory Elements 


After the archaic and AMH lineages diverged, each accu- 
mulated genetic variation specific to each group. Variants 
fixed in each lineage are likely to be enriched in genomic re- 
gions that influence traits that experienced positive selec- 
tion. We tested whether human- and archaic-specific 
fixed variants are enriched compared with other variants 
in circadian genes and their promoters and in annotated 
cCREs within 1 Mb (fig. 2). We found that human- and 
archaic-specific fixed variants are enriched in circadian 
genes (Fisher's exact test; human: odds ratio [OR] = 1.84, 
P=7.06e—12; archaic: OR=1.13, P=0.023) and distal 
regulatory elements (Fisher's exact test; human: OR= 
1.25, P=8.39e—4: archaic: OR = 1.16, P= 6.15e—5) com- 
pared with variants derived on each lineage but not fixed. 
Promoter regions have a similar enrichment pattern as 
that in gene and regulatory regions, but the P-values are 
high (Fisher's exact test; human: OR = 1.21, P= 0.65; ar- 
chaic: OR = 1.09, P= 0.63). This is likely due to the small 
number of such variants in promoters. These results sug- 
gest that both groups had a greater functional divergence 
in genomic regions related to circadian biology than 
expected. 


Several Core Circadian Genes Have Evidence of 
Alternative Splicing between Humans and 
Archaic Hominins 


We find only two archaic-specific coding variants in circa- 
dian genes: one missense and one synonymous. The mis- 
sense variant (hg19: chr17_46923411_A_G) is in the 
gene CALCOCOZ, calcium-binding and coiled-coil domain- 
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Fic. 1.—Did the sharing of functionally diverged alleles from archaic hominins influence human circadian biology? (A) AMHs and archaic hominins 
evolved separately at different latitudes for hundreds of thousands of years. The ancestors of modern Eurasian humans left Africa approximately 70 ka 
and admixed with archaics, likely in southwestern Asia. The shaded range in Eurasia represents the approximate Neanderthal range. The dot in Asia represents 
the location of the sequenced Denisovan individual in the Altai Mountains; the full range of Denisovans is currently unknown. Silhouettes from phylopic.org. 
(B) After the split between the human and archaic lineages, each group accumulated variation and evolved in their respective environments for approximately 
700 ka. We first test for evidence for divergent circadian evolution during this time. Humans acquired introgressed alleles from Neanderthals and from 
Denisovans around 60 and 45 ka, respectively. These alleles experienced strong selective pressures; however, ~40% of the genome retains archaic ancestry 
in some modern populations. The second question we explore is whether introgression made contributions to human circadian biology. (C) The core circadian 
clock machinery is composed of several transcription factors (ovals) that dimerize and interact with E-box enhancer elements and each other to create a nega- 
tive feedback loop. We defined a set of 246 circadian genes through a combination of literature search, expert knowledge, and existing annotations 
(supplementary table S1 and fig. S1, Supplementary Material online; Materials and Methods). Lines with arrows represent activation, and lines with bars re- 
present suppression. 


containing protein 2. SIFT, PolyPhen, and CADD all predict analysis (the Altai, the Vindija, the Chagyrskaya 


that the variant does not have damaging effects. The se- 
cond variant (hg19: chr7_119914770_G_T) is in the gene 
KCND2, which encodes a component of a voltage-gated 
potassium channel that contributes to the regulation of 
the circadian rhythm of action potential firing, but it is syn- 
onymous and the variant effect predictors suggest it is 
tolerated. 

To explore potential splicing differences in circadian 
genes between humans and archaics, we applied SpliceAl 
to predict whether any sequence differences between 
modern humans and archaics are likely to modify splicing 
patterns. Four archaic individuals were included in this 


Neanderthals, and the Altai Denisovan) (Meyer et al. 
2012; Prüfer et al. 2014, 2017; Mafessoni et al. 2020). 
We found that 28 genes contained at least one archaic- 
specific variant predicted to result in alternative splicing in 
archaics. These included several of the core clock genes 
CLOCK, PER2, RORB, RORC, and FBXL13 (fig. 3A and G; 
supplementary table S3, Supplementary Material online). 
For example, the variant (hg19: chr2:239187088- 
2391870839) in the first intron of PER2 is predicted to result 
in a longer 5’ UTR. The splice-altering variants were largely 
specific to the two different archaic lineages (fig. 3A), with 
13 specific to the Denisovan, eight shared among the three 
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Enrichment 
Human specific 1.25 (247) 1.21 (6) 1.84 (156) 
Archaic specific 1.16 (807) 1.09 (19) 1.13 (341) 


Regulatory Promoters Genes 
Elements 


Circadian genes 


Fic. 2.—Human- and archaic-specific fixed variants are enriched in cir- 
cadian regulatory, promoter, and gene regions. Human-specific fixed var- 
iants are significantly enriched compared with variants that are not fixed in 
circadian regulatory elements (Fisher's exact: OR = 1.25, P= 8.39e—4) and 
gene regions (Fisher's exact: OR = 1.84, P= 7.06e—12). Promoters show a 
similar enrichment, but the higher P-value is the result of the small number 
of variants (Fisher's exact test: OR = 1.21, P=0.65). Likewise, archaic- 
specific variants are enriched in circadian regulatory regions (Fisher's exact: 
OR = 1.16, P=6.15e—5) and gene regions (Fisher's exact: OR = 1.13, P= 
0.023), with the promoters showing a similar trend (Fisher's exact test: OR 
= 1.09, P= 0.63). The numbers in parentheses give the counts of fixed var- 
iants observed in each type of element. Regulatory elements were defined 
based on the ENCODE cCREs. 


Neanderthals, and only one shared among all four archaic 
individuals. 


Circadian Gene Regulatory Divergence between 
Humans and Archaic Hominins 


Given the enrichment of variants in regulatory regions of 
circadian genes, we sought to explore the potential for dif- 
ferences in circadian gene regulation between humans and 
archaics with causes beyond single lineage-specific var- 
iants. We leveraged an approach we recently developed 
for predicting gene regulatory differences between modern 
and archaic individuals from combinations of genetic var- 
iants (Colbran et al. 2019). The approach uses PredixXcan, 
an elastic net regression method, to impute gene transcript 
levels in specific tissues from genetic variation. Previous 
work demonstrated that this approach has a modest de- 
crease in performance when applied to Neanderthals, but 
that it can accurately applied between humans and 
Neanderthals for thousands of genes. Here, we quantify 
differences in the predicted regulation of the 246 circadian 
genes between 2,504 humans in the 1000 Genomes 
Project (1000 Genomes Project Consortium 2010) and 
the archaic hominins. The predicted regulation values are 
normalized to the distribution in the training set from the 
Genotype Tissue Expression Atlas (GTEx). 

We first analyzed gene regulation predictions in the core 
circadian clock genes. Archaic gene regulation was at the 
extremes of the human distribution for many core clock 
genes including PER2, CRY1, NPAS2, RORA, and NR1D1 


(fig. 4; supplementary fig. S2, Supplementary Material on- 
line). For example, the regulation of PER2 in the Altai and 
Vindija Neanderthals is lower than 2,491 of the 2,504 
(99.48%) modern humans considered. The Denisovan has 
a predicted PER2 regulation that is lower than 2,410 
(96.25%). Expanding to all circadian genes and requiring ar- 
chaic regulation to be more extreme than all humans 
(Materials and Methods), we identified 24 circadian genes 
across 23 tissues with strong divergent regulation between 
humans and at least one archaic hominin (fig. 38; 
supplementary table S4, Supplementary Material online). For 
example, all archaic regulation values for RORA, a core clock 
gene, are lower than for any of the 2,504 modern humans. 
We found that 16 of these genes (supplementary fig. S3 
and table S4, Supplementary Material online), including 
RORA, MYBBP1A, and TIMELESS, were divergently regulated 
(DR) in all archaic individuals. This represents 6.5% of all the 
circadian genes (fig. 3C). Surprisingly, the two Neanderthals 
only shared one DR gene not found in the Denisovan, while 
the Altai Neanderthal and Denisovan shared seven not found 
in Vindija (fig. 3B). The Altai and Vindija Neanderthals re- 
present deeply diverging lineages, and this result suggests 
that they may have experienced different patterns of diver- 
gence in the regulation of their circadian genes. 

Given these differences in circadian gene regulation be- 
tween humans and archaics, we tested whether circadian 
genes are more likely to be DR than other gene sets. Each ar- 
chaic individual shows nominal enrichment for divergent regu- 
lation of circadian genes, and the enrichment was stronger 
(~1.2x) in the Altai Neanderthal and Denisovan individuals. 
However, given the small sample size, the P-values are moder- 
ate (permutation test; Altai: OR = 1.21, P=0.19; Vindija: OR 
= 1.05, P=0.43; Denisovan: OR= 1.20, P=0.24). 


Did Introgressed Archaic Variants Influence Modern 
Human Circadian Biology? 


The previous sections demonstrate lineage-specific genetic 
variation in many genes and regulatory elements essential 
to the function of the core circadian clock and related path- 
ways. Given this evidence of functional differences be- 
tween archaic hominins and AMH in these systems, we 
next evaluated the influence of archaic introgression on 
AMH circadian biology. 


Introgressed Variants Are Enriched in Circadian 
Gene eQTL 


Given the differences between archaic and modern se- 
quences of circadian genes and their regulatory elements, 
we investigated whether Neanderthal introgression con- 
tributed functional circadian variants to modern Eurasian 
populations. We considered a set of 863,539 variants 
with evidence of being introgressed from archaic hominins 
to AMH (Browning et al. 2018). These variants were 
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Fic. 3.—Many circadian genes have evidence of alternative splicing and divergent regulation between modern and archaic hominins. (A) The distribution 
of the 28 predicted archaic-specific splice-altering variants (SAVs) in circadian genes across archaic individuals. Most are specific to either the Denisovan or 
Neanderthal lineage (supplementary table $3, Supplementary Material online). (8) The sharing of predicted DR gene/tissue pairs across three archaic indivi- 
duals. Predictions were not available for the Chagyrskaya Neanderthal. Seventeen DR gene/tissue pairs were present in all three archaics (representing 16 
unique genes). Additionally, seven gene/tissue DR pairs are shared between the Altai Neanderthal and the Denisovan individual. One pair is shared between 
the Vindija Neanderthal and the Denisovan (supplementary table S4, Supplementary Material online). (C) The proportion of circadian genes containing archaic 
SAVs predicted by SpliceAl (11.4%) or DR circadian genes predicted by PredixXcan (6.5%). Thus, 17.9% of the circadian genes are predicted to contain differ- 


ences to AMH via these mechanisms. 


identified using the Sprime algorithm, which searches for 
regions containing a high density of alleles in common with 
Neanderthals and not present or at very low frequency in 
Africans. Since many approaches have been developed to 
identify introgressed variants, we also considered two stricter 
sets: 47,055 variants that were supported by all of six different 
introgression maps (Sankararaman et al. 2014; Vernot et al. 
2016; Browning et al. 2018; Steinrücken et al. 2018; Skov 
et al. 2020; Schaefer et al. 2021) and 755,653 variants that 
were supported by Sprime and at least one other introgression 
map. As described below, our main results replicated on both 
of these stricter sets. 

We first tested whether the presence of introgressed var- 
iants across modern individuals is associated with the ex- 
pression levels of any circadian genes, that is, whether 
the introgressed variants are expression quantitative trait 
loci (eQTL). We identified 3,857 introgressed variants asso- 
ciated with the regulation of circadian genes in modern 
non-Africans (supplementary table S5, Supplementary 
Material online). The genes PTPRJ, HTR1B, NR1D2, 
CLOCK, and ATOH7 had the most eQTL (304, 273, 262, 
256, and 252, respectively). We found introgressed circa- 
dian eQTL for genes expressed in all tissues in GTEx, except 
the kidney cortex. Notably, several of these circadian genes 
(e.g., NR1D2 and CLOCK) with introgressed eQTL were also 
found to be DR in our comparison of modern and archaic 
gene regulation. This indicates that some of the archaic- 
derived variants that contributed to divergent regulation 
were retained after introgression and continue to influence 
circadian regulation in modern humans. 


Introgressed variants are significantly more likely to be 
eQTL for circadian genes than expected by chance from 
comparison to all eQTL (fig. 5A; Fisher's exact test: OR = 
1.45, P=9.71e—101). The stricter set of introgressed var- 
iants identified by Browning et al. plus at least one other 
introgression map had similar levels of eQTL enrichment 
for circadian genes (OR = 1.47; P= 2.4e—103). The highest 
confidence set of introgressed variants that were identified 
by all six maps considered had even stronger enrichment 
(OR = 1.68; P=6.5e—23). 

Most core circadian genes are expressed broadly across 
tissues; the fraction expressed in each GTEx tissue ranges 
from 57% (whole blood) to 83% in testis and an average 
of 72% (supplementary table S6, Supplementary Material 
online). As a result, we anticipated that the enrichment of in- 
trogressed variants among eQTLs for circadian genes would 
hold across tissues. Examining the associations in each tissue, 
we found that introgressed eQTL showed significant enrich- 
ment for circadian genes in most tissues (34 of 49; fig. 5B; 
supplementary table S7, Supplementary Material online) 
and trended this way in all but five. Given that tissues in 
GTEx have substantial differences in sample size and cellular 
heterogeneity, statistical power to detect enrichment differs. 
We anticipate that this is the main driver of differences in en- 
richment across tissues. 

These results suggest that circadian pressures were wide- 
spread across tissues. Given the previously observed depletion 
for introgressed variants in regulatory elements and eQTL 
(Petr et al. 2019; Rinker et al. 2020; Telis et al. 2020), this en- 
richment for circadian genes among introgressed eQTL is 
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Fic. 4.—Many circadian genes are divergently regulated between modern humans and archaic hominins. Comparison of the imputed regulation of core 
circadian genes between 2,504 humans in 1000 Genomes Phase 3 (histogram) and three archaic individuals (vertical lines). For each core circadian gene, the 
tissue with the lowest average P-value for archaic difference from humans is plotted. Archaic gene regulation is at the extremes of the human distribution for 
several core genes: CRY7, PER2, NPAS2, NR1D1, and RORA. See supplementary figure S2, Supplementary Material online for all core clock genes and 
supplementary figure S3, Supplementary Material online for all DR circadian genes. 


surprising and suggests that the archaic circadian alleles 
could have been beneficial after introgression. 


Introgressed Variants Predominantly Increase Propensity 
for Morningness 


After observing that circadian gene expression is influenced 
by archaic variants, we evaluated whether these effects are 
likely to result in a change in organism-level phenotype. To 
do this, we evaluated evidence that introgressed variants 
influence chronotype. The heritability of chronotype has 
been estimated in a range from 12% to 38% (Jones et al. 
2016, 2019; Lane et al. 2016), and previous studies have 
identified two introgressed loci associated with sleep pat- 
terns (Dannemann and Kelso 2017; Putilov et al. 2019). 
We recently found modest enrichment for heritability of 
chronotype (morning/evening person phenotype in a 
GWAS of the UK Biobank) among introgressed variants 
genome wide using stratified linkage disequilibrium (LD) 
score regression (heritability enrichment: 1.58, P=0.25) 
(McArthur et al. 2021). This analysis also suggested that intro- 
gressed variants were more likely to increase morningness. 


To test for this proposed directional effect, we calculated 
the cumulative fraction of introgressed loci associated with 
chronotype in the UK Biobank that increase morningness 
(after collapsing based on LD at R*>0.5 in EUR). The 
introgressed loci most strongly associated with chronotype 
increase the propensity for morningness (fig. 6; 
supplementary tables S8 and S9, Supplementary Material 
online). As the strength of the association with morning- 
ness decreases, the bias begins to decrease, but the effect 
is maintained well past the genome-wide significance 
threshold (P < 5e—8). When focusing the analysis on intro- 
gressed variants in proximity (<1 Mb) to circadian genes, 
the pattern becomes even stronger. The bias toward morn- 
ingness remains above 80% at the genome-wide signifi- 
cance threshold. This result also held when limiting to 
introgressed variants found in Browning plus one or all 
other introgression maps considered (supplementary fig. 
S4, Supplementary Material online). This suggests that in- 
trogressed variants act in a consistent direction on chrono- 
type, especially when they influence circadian genes. 

Circadian rhythms are involved in a wide variety of bio- 
logical systems. To explore other phenotypes potentially 
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influenced by the introgressed circadian variants, we 
evaluated evidence for pleiotropic associations. First, we re- 
trieved all the genome-wide associations reported for intro- 
gressed variants in the Open Targets Genetics (https:/ 
genetics.opentargets.org) database, which combines 
GWAS data from the GWAS Catalog, UK Biobank, and sev- 
eral other sources. Introgressed circadian variants are asso- 
ciated with traits from a diverse range of categories 
(supplementary table $10, Supplementary Material online). 
Associations with blood-related traits are by far the most 
common; however, this is likely because they have more 
power in the UK Biobank. Overall, circadian introgressed 
variants are significantly more likely to have at least one 
trait association than introgressed variants not in the circa- 
dian set (Fisher's exact test: OR=1.25, P=7.03e—25) 
(supplementary fig. SSA, Supplementary Material online). 
The circadian variants also associate with significantly 
more traits per variant than the noncircadian set (Mann- 
Whitney U test: P= 9.93e—14) (supplementary fig. S58 
and table S11, Supplementary Material online). These 
results suggest effects for introgressed circadian variants 
beyond chronotype. 


Evidence for Adaptive Introgression at Circadian Loci 


The gene flow from Eurasian archaic hominins into AMH 
contributed to adaptations to some of the new environ- 
mental conditions encountered outside of Africa (Racimo 
et al. 2015). The above analyses demonstrate the effects 
of introgressed variants on circadian gene regulation and 
chronotype. To explore whether these circadian regions 
show evidence of adaptive introgression, we considered 
three sets of introgressed regions predicted to have contrib- 
uted to AMH adaptation: one from an outlier approach 
based on allele frequency statistics (Racimo, Marnetto, 
and Huerta-Sanchez 2017) and two from recent machine 
learning algorithms: genomatnn (Gower et al. 2021) and 
MalAdapt (Zhang et al. 2023). We intersected the circadian 
introgressed variants with the adaptive introgression re- 
gions from each method. 

We identified 47 circadian genes with evidence of adap- 
tive introgression at a nearby variant from at least one of 
the methods (supplementary table S12, Supplementary 
Material online). No region was supported by all three 
methods; however, six were shared between Racimo 
and MalAdapt and three were shared by Racimo and 
genomatnn. The relatively small overlap between these 
sets underscores the challenges of identifying adaptive 
introgression. Nonetheless, these represent promising can- 
didate regions for further exploration of the effects of intro- 
gressed variants on specific aspects of circadian biology. For 
example, an introgressed haplotype on chr10 tagged by 
rs76647913 was identified by both MalAdapt and 
Racimo. This introgressed haplotype is an eQTL for the 
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nearby ATOH7 gene in many GTEx tissues. ATOH7 is a cir- 
cadian gene that is involved in retinal ganglion cell develop- 
ment, and mice with this gene knocked out are unable to 
entrain their circadian clock based on light stimuli 
(Brzezinski et al. 2005). 


Latitudinal Clines for Introgressed Circadian Loci 


Motivated by the previous discovery of an introgressed 
haplotype on chr2 that is associated with chronotype and 
increases in frequency with latitude (Dannemann and 
Kelso 2017; Putilov et al. 2019), we also tested each intro- 
gressed circadian variant for a correlation between allele 
frequency and latitude in modern non-African populations 
from the 1000 Genomes Project. 

The strongest association between latitude and fre- 
quency was a large chromosome 2 haplotype that contains 
the previously discovered introgressed single nucleotide 
polymorphism (SNP) (rs75804782, R=0.85) associated 
with the chronotype. This haplotype is present in all 
non-African populations, and rs61332075 showed the 
strongest latitudinal cline (R = 0.87). The second strongest 
consisted of a smaller haplotype of introgressed variants a 
few kilobases of upstream of the previous haplotype 
(tagged by rs35333999 and rs960783) that overlaps the 
core circadian gene PER2. These variants have a correlation 
between latitude and frequency of ~0.68. They are also in 
moderate LD (R? of ~0.35 in EUR) with an additional intro- 
gressed variant (rs62 194932) that has a similar latitudinal 
cline of 0.70 (supplementary fig. S6 and table S13, 
Supplementary Material online). These variants are each 
in very low LD with the previously discovered haplotype 
(R? of ~0.01) and are each supported by multiple introgres- 
sion maps. Moreover, these introgressed variants are 
absent in all EAS populations, absent or at very low fre- 
quency in SAS (<3%), and at higher frequency in EUR po- 
pulations (~13%). 

The EUR-specific introgressed variant rs35333999 causes a 
missense change in the PER2 protein (V903I) that overlaps a 
predicted interaction interface with PPARG. PER2 controls li- 
pid metabolism by directly repressing PPARG’s proadipogenic 
activity (Grimaldi et al. 2010). The rs62194932 variant is an 
eQTL of HES6 in the blood in the eQTLGen cohort (Võsa et 
al. 2021). HES6 encodes a protein that contributes to circa- 
dian regulation of LDLR and cholesterol homeostasis (Lee et 
al. 2012). 

Thus, this genomic region that includes circadian genes 
and introgressed variants associated with chronotype has 
a population-specific structure and at least two distinct 
sets of introgressed variants with latitudinal clines and func- 
tional links to lipid metabolism. PER2 is also predicted to 
have lower gene regulation in archaic hominins than 
most humans (fig. 4), and the Vindija Neanderthal carries 
a lineage-specific variant in this gene that has 
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splice-altering effects. These results together suggest that 
PER2 may have experienced multiple functional changes 
in different modern and archaic lineages, with potential 
adaptive effects mediated by introgression. 

We did not discover any other significant associations 
between latitude and frequency for other introgressed cir- 
cadian loci. The rapid migration and geographic turnover of 
populations in recent human history are likely to obscure 
many latitude-dependent evolutionary signatures, so we 
did not anticipate many circadian loci would have a strong 
signal. 


Discussion 


The Eurasian environments where Neanderthals and 
Denisovans lived for several hundred thousand years are lo- 
cated at higher latitudes with more variable photoperiods 
than the landscape where AMH evolved before leaving 
Africa. Evaluating genetic variation that arose separately 
in each of the archaic and AMH lineages after their split 
~700 Ma, we identified lineage-specific genetic variation 
in circadian genes, their promoters, and flanking distal 
regulatory elements. We found that both archaic- and 
human-specific variants are observed more often than ex- 
pected in each class of functional region. This result sug- 
gests that, while each group evolved separately during 
hundreds of thousands of years in divergent environments, 
both experienced pressure on circadian-related variation. 
Leveraging sequence-based machine learning methods, 
we identified many archaic-specific variants likely to influ- 
ence circadian gene splicing and regulation. For example, 
core clock genes (CLOCK, PER2, RORB, RORC, and FBXL13) 
have archaic variants predicted to cause alternative splicing 
compared with AMH. Several core genes were also predicted 
in archaics to be at the extremes of human gene regulation, 
including PER2, CRY1, NPAS2, RORA, and NR1D1. 
Surprisingly, the Altai Neanderthal shared more divergent 
regulation in the circadian genes with the Denisovan individ- 
ual than the Vindija Neanderthals. The two Neanderthals re- 
present populations that were quite distantly diverged with 
substantially different histories and geographical ranges. The 
Denisovan and Altai Neanderthal also come from the same re- 
gion in Siberia, while the Vindija Neanderthal came from a re- 
gion in Croatia with slightly lower latitude. 

Introgression introduced variation that first appeared in 
the archaic hominin lineage into Eurasian AMH. While 
most of this genetic variation experienced strong negative 
selection in AMH, a smaller portion is thought to have pro- 
vided adaptive benefits in the new environments (Racimo 
et al. 2015). Given the divergence in many circadian genes’ 
regulation, we explored the landscape of introgression on 
circadian genes. We first looked at introgressed circadian 
variants that are likely to influence gene regulation in 
AMH. Variants in this set are observed more often than 


expected, suggesting the importance of maintaining circa- 
dian variation in the population. We also verified that these 
results held over variants identified by different methods for 
calling archaic introgression. 

We then evaluated the association of these introgressed 
variants with variation in the circadian phenotypes of 
Eurasians. We previously reported a modest enrichment 
among introgressed variants for heritability of the morn- 
ing/evening person phenotype (McArthur et al. 2021). 
Here, we further discovered a consistent directional effect 
of the introgressed circadian variants on chronotype. The 
strongest associated variants increase the probability of 
being a morning person in Eurasians. 

While it is not immediately clear why increased morning- 
ness would be beneficial at higher latitudes, considering 
this directional effect in the context of clock gene regula- 
tion and the challenge of adaptation to higher latitudes 
suggests an answer. In present-day humans, behavioral 
morningness is correlated with a shortened period of the 
circadian molecular clockworks in individuals. This earlier 
alignment of sleep/wake with external timing cues is a con- 
sequence of a quickened pace of the circadian gene net- 
work (Brown et al. 2008). Therefore, the morningness 
directionality of introgressed circadian variants may indi- 
cate selection toward a shortened circadian period in the 
archaic populations living at high latitudes. Supporting 
this interpretation, shortened circadian periods are required 
for synchronization to the extended summer photoperiods of 
high latitudes in Drosophila, and selection for shorter periods 
has resulted in latitudinal clines of decreasing period with in- 
creasing latitude, as well as earlier alignment of behavioral 
rhythms (Hut et al. 2013). In addition, Drosophila populations 
exhibit decreased amplitude of behavioral rhythms at higher 
latitudes, which is also thought to aid in synchronization to 
long photoperiods (Hut et al. 2013). 

Our finding that introgressed circadian variants generally 
decrease gene regulation of circadian genes suggests that 
they could lead to lower amplitude clock gene oscillations. 
However, when assayed in present-day humans, there is 
not a strong correlation between the overall expression le- 
vel of NR1D1 and the transcriptional amplitudes of other 
clock genes within individuals (Brown et al. 2008), and 
quantitative modeling of the mammalian circadian clock- 
works suggests that stable clock gene rhythms can result 
across a wide range of absolute levels of gene expression 
as long as the stoichiometric ratios of key positive and nega- 
tive clock genes are reasonably conserved (Kim and Forger 
2012). Interestingly, the lower transcriptional amplitude of 
NR1D1 does confer greater sensitivity of the present-day 
human clockworks to resetting stimuli, a potentially adap- 
tive characteristic for high latitudes (Brown et al. 2008). 

Thus, given the studies of latitudinal clines and adapta- 
tion from Drosophila and the nascent understanding of 
clock gene contributions to behavioral phenotypes in 
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present-day humans, the directional effects of introgressed 
circadian gene variants toward early chronotype and de- 
creased gene regulation we observed can be viewed as po- 
tentially adaptive. More complex chronotype phenotyping 
and mechanistic studies of the variants of interest are 
needed to fully understand these observations. 

Finally, to explore evidence for positive selection on in- 
trogressed variants in AMH, we analyzed results from three 
recent methods for detecting adaptive introgression. All 
methods identified circadian loci as candidates for adaptive 
introgression. However, we note that the predictions of 
these methods have only modest overlap with one another, 
underscoring the difficulty of identifying adaptive intro- 
gression. Nonetheless, many of these loci, especially those 
supported by both Racimo and MalAdapt, are good candi- 
dates for adaptive introgression given their functional asso- 
ciations with circadian genes 

Several limitations must be considered when interpret- 
ing our results. First, it is challenging to quantify the com- 
plexity of traits with a large behavioral component (like 
chronotype) and infer their variation from genomic infor- 
mation alone. Nevertheless, we believe our approach of fo- 
cusing on molecular aspects (splicing and gene regulation) 
of genomic loci with relevance to circadian biology, in par- 
allel to GWAS-based associations, lends additional support 
to the divergence in chronotype between archaic hominins 
and modern humans. Second, we also note that circadian 
rhythms contribute to many biological systems, so the var- 
iants in these genes tend to be associated with a variety of 
phenotypes. Thus, there is also the potential that selection 
acted on other phenotypes influenced by circadian vari- 
ation than those related directly to chronotype. Third, given 
the complexity of circadian biology, there is no gold stand- 
ard set of circadian genes. We focus on the core clock genes 
and a broader set of expert-curated genes relevant to circa- 
dian systems, but it is certainly possible that other genes 
with circadian effects are not considered. Fourth, recent 
adaptive evolution is challenging to identify, and this is es- 
pecially challenging for introgressed loci. Nonetheless, we 
find several circadian loci with evidence of adaptive intro- 
gression from more than one method. Finally, given the 
many environmental factors that differed between 
African and non-African environments, it is difficult to de- 
finitively determine whether selection on a particular locus 
was the result of variation in light levels versus other related 
factors, such as temperature. Nonetheless, given the ob- 
served modern associations with chronotype for many of 
these variants, we believe it is a plausible target. 

In conclusion, studying how humans evolved in the face 
of changing environmental pressures is necessary to under- 
standing variation in present-day phenotypes and the 
potential trade-offs that influence the propensity to differ- 
ent diseases in modern environments (Benton et al. 2021). 
Here, we show that genomic regions involved in circadian 
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biology exhibited substantial functional divergence be- 
tween separate hominin populations. Furthermore, we 
show that introgressed variants contribute to variation in 
AMH circadian phenotypes today in ways that are consist- 
ent with an adaptive benefit. 


Materials and Methods 


Circadian Gene Selection 


Circadian biology is a complex system due to its high im- 
portance in the functioning of biological timing in diverse 
biological systems. For that reason, determining which 
genes are crucial for selection to environmental response 
related to light exposure is not a straightforward process. 
To address this issue, we look at different sources of gen- 
ome annotation databases and searched for genes and var- 
iants associated with circadian-related phenotypes. We 
considered all human protein-coding genes in the Gene 
Ontology database annotated with the GO: 0007623 
(“circadian rhythm”) term or terms annotated with the re- 
lationship “is_a,” “part_of,” “occurs_in,” or “regulates” 
circadian rhythm. We also considered genes containing ex- 
perimental or orthologous evidence of circadian function in 
the Circadian Gene Database (CGDB), the GWAS Catalog 
genes containing “chronotype” or “circadian rhythm” as- 
sociated variants, and a curated set of genes available in 
WikiPathways (Martens et al. 2021). The final set of circa- 
dian genes was curated by Dr. Douglas McMahon. 

To select the candidate circadian genes with the highest 
confidence, we defined a hierarchy system where genes an- 
notated by McMahon or annotated in three out of four 
other sources receive a “high” level of confidence. Genes 
with evidence from two out of four of the sources are as- 
signed a “medium” level of confidence. Genes annotated 
as circadian only in one out of four sources are assigned 
to “low” level of confidence and not considered in our cir- 
cadian gene set. We then defined our set of circadian var- 
iants from the 1000 Genomes Project using the list of 
circadian genes. The variants are included in the analysis 
of coding, noncoding, regulatory, eQTL, human-specific, 
archaic-specific, and introgressed variants. 


wou 


Definition of Lineage-Specific Variants 


To identify candidate variants that are specific to the human 
and the archaic lineages, we used a set of variants pub- 
lished by Kuhlwilm and Boeckx (2019). The variants were 
extracted from the high-coverage genomes of three archa- 
ics: a 122,000-year-old Neanderthal from the Altai Mountains 
(52x coverage), a 52,000-year-old Neanderthal from Vindija in 
Croatia (30x coverage), and a 72,000-year-old Denisovan 
from the Altai Mountains (30x coverage). The variants were 
called in the context of the human genome hg19/GRCh37 ref- 
erence. The total number of variant sites after applying filters 
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for high coverage sites and genotype quality is 4,437,803. A 
human-specific variant is defined as a position where all the 
humans in the 1000 Genomes Project carry the derived allele 
and all the archaics carry the ancestral allele. An archaic specific 
is defined as a position where all the archaics carry the derived 
allele and the derived allele is absent or extremely rare 
(<0.00001) across all human populations. Note that intro- 
gressed archaic alleles are not included in the “archaic- 
specific” set. These criteria resulted in 9,424 human-specific 
and 33,184 archaic-specific variants. 


Enrichment of Lineage-Specific Variants among 
Functional Regions of the Genome 


We intersected the sets of lineage-specific variants with sev- 
eral sets of annotated functional genomic regions. Inside cir- 
cadian gene regions (Gencode v29), we found 156 
human-specific variants and 341 archaic-specific variants. In 
circadian promoter regions, we found six human-specific var- 
iants and 19 archaic-specific variants. Promoters were defined 
as regions 5-kb upstream to 1-kb downstream from a tran- 
scription start site (TSS). In distal regulatory elements, we 
found 247 human-specific variants and 807 archaic-specific 
variants. For this last set, we considered cCREs published by 
ENCODE (Moore et al. 2020) within 1 Mb of the circadian 
genes. To compute whether lineage-specific variants are 
more abundant than expected in circadian genes, we applied 
a Fisher's exact test to the sets of human- and archaic-specific 
variants in regulatory, promoter, and gene regions. 


Genes Containing Archaic Variants with Evidence of 
Alternative Splicing 


We used a set of archaic variants annotated with the 
splice-altering probabilities to identify circadian genes 
that may be differentially spliced between archaic hominins 
and AMH (Brand et al. 2023). We considered variants from 
four archaic individuals: the Altai, Chagyrskaya, and Vindija 
Neanderthals and the Altai Denisovan. These archaic var- 
iants were annotated using SpliceAl (Jaganathan et al. 
2019), and we considered any variant with a maximum del- 
ta or splice-altering probability > 0.2. We identified 36 
archaic-specific splice-altering variants, defined as those 
variants absent from 1000 Genomes Project, among 28 
circadian genes. Next, we tested for enrichment among 
this gene set using an empirical null approach (McArthur 
t al. 2022; Brand et al. 2023). We shuffled the maximum 
eltas among 1,607,350 variants 10,000 times and 
counted the number of circadian genes with a splice- 
altering variant in each iteration. Enrichment was calculated 
as the number of observed genes (N = 28) divided by the 
mean gene count among 10,000 shuffles. In addition to 
all genes with archaic-specific variants, we considered six 
other subsets among these variants: 1) genes with variants 
private to the Altai Neanderthal, 2) genes with variants 
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private to the Chagyrskaya Neanderthal, 3) genes with var- 
iants private to the Altai Denisovan, 4) genes with variants 
private to all Neanderthals, 5) genes with variants shared 
among all archaic individuals, and 6) genes with variants 
private to the Vindija Neanderthal. Finally, we considered 
a subset of splice-altering variants that were identified as 
tag SNPs by Vernot et al. (2016). 


Predixcan 


To understand the difference in circadian biology between 
present-day humans and archaic hominins, we analyzed 
predictions on gene regulation. We considered the results 
from PrediXcan gene regulation predictions across 44 tis- 
sues from the PredictDB Data Repository (http:/predictdb. 
org/). The models were trained on GTEx V6 using variants 
identified in 2,504 present-day humans in the 1000 
Genomes Project Phase 3 within 1 Mb of each circadian 
gene. The original analysis includes predictions for 17,748 
genes for which the models explained a significant amount 
of variance in gene expression in each tissue (false discovery 
rate < 0.05). The prediction models were also applied to the 
Altai and Vindija Neanderthals and the Denisovan. The result- 
ing predictions are normalized values of the distribution ob- 
served in GTEx individuals used to train the original 
prediction models. Each prediction contains an empirical 
P-value, which was calculated for each gene and tissue pair 
to define genes that are DR between archaic hominins and 
humans. The P-value is obtained by calculating the proportion 
of humans from the 1000 Genomes Project that have predic- 
tions more extreme compared to the human median than the 
archaic individual. Significantly, DR genes are defined as those 
where the archaic prediction falls outside the distribution of 
humans in the 1000 Genomes Project predictions. 

We tested whether the circadian genes in our set are 
more likely to be DR compared with an empirical null distri- 
bution from random gene sets of the same size. We ac- 
count for the fact that some genes are modeled in more 
tissues than others by matching the distribution of tissues in 
which each gene could be modeled in the random sets to 
our set. Among 1,467 DR genes in the Altai Neanderthal, 
we find 23 DR circadian genes out of the total 236 genes in 
the circadian set. We iterate through the permutation analysis 
1,000,000 times and find an enrichment of 1.21 (P=0.19). A 
similar analysis is done in the Vindija Neanderthal (1,536 tota 
DR, 21 circadian DR, enrichment of 1.05, P=0.43) and the 
Denisovan individual (1,214 total DR, 19 circadian DR, enrich- 
ment of 1.20, P=0.24). In this study, we define a set of DR 
genes as the intersection between DR genes in all three archa- 
ics, resulting in a set of 16 genes. 


Enrichment of Introgressed Variants in eQTL 


We performed an enrichment analysis using Pearson's 
chi-squared test to evaluate if there is over-representation 
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of introgressed alleles in our set of circadian variants using 
the GTEx dataset. We did a liftOver of the GTEx v8 data set 
from hg38 to hg19. The original hg38 set contains 
4,631,659 eQTLs across 49 tissues. After the LiftOver, 
4,608,446 eQTLs remained, with the rest not mapping. 
We used the archaic introgressed variants data set from 
Browning et al. (2018). The set contains 863,539 variants 
that are introgressed in humans originating in archaic ho- 
minins. We performed an intersection between the set of 
genes containing evidence for eQTLs and our set of 246 cir- 
cadian genes to retrieve a subset of variant sites with evi- 
dence of being eQTL in circadian genes. The resulting 
subset contained 97,441 circadian eQTLs in 49 tissues 
and 239 genes. We further intersected the introgressed 
variants and the set of eQTL, resulting in 128,138 intro- 
gressed eQTLs. The final set of eQTLs that are circadian 
and also introgressed is 3,857. 


Direction of Effect of Chronotype Associations 


To explore the effect of archaic introgression in circadian 
genes on human chronotype, we quantified the direction 
of the effect of variants associated to a morning/evening 
person trait in a GWAS analysis of the UK Biobank (http:/ 
www.nealelab.is/uk-biobank/). The variants were LD 
clumped using PLINK v1.9 (R? > 0.5). We generated cumu- 
lative proportion values on the beta values assigned to each 
associated variant on an ascending order of P-values. 


Detection of Latitudinal Clines in Chronotype 
Associations 


To evaluate latitudinal clines in chronotype-associated var- 
iants, we assigned a latitude to each of the Eurasian 1000 
Genomes Project populations. The latitude of diaspora 
populations was set to their ancestral country (GIH 
Gandhinagar in Gujarat: 23.223; STU Sri Jayawardenepura 
Kotte: 6.916667; ITU Amaravati in Andhra Pradesh: 
16.5131; CEU: 52.372778). CEU was assigned a latitude in 
Amsterdam, following an analysis that shows that this group 
is more closely related to Dutch individuals (Lao et al. 2008). 
We then used the LDlink API to retrieve allele frequencies for 
each introgressed morningness variant in Eurasian indivi- 
duals (Machiela and Chanock 2015). Variants that follow 
a latitudinal cline were identified using linear regression sta- 
tistics requiring correlation coefficient (R > 0.65) and P-value 
(P<0.5). 


Detection of Pleiotropy in the Set of Introgressed 
Circadian Variants 


To understand the extent of different phenotypes asso- 
ciated with the introgressed circadian variants, we first ex- 
tracted genome-wide associations from Open Targets 
Genetics (https:/genetics.opentargets.org/) for each of 
the variants with evidence of introgression (Browning et 
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al. 2018). Only the variants with significant P-values were 
analyzed. The P-value threshold was set at the genome- 
wide significance level (P=5e—8). We split the variants 
into two sets: introgressed circadian and introgressed non- 
circadian. Many of these variants are not associated with 
any phenotype. We performed a Fisher's exact test to ana- 
lyze which of the two sets contains a higher ratio of SNPs 
with at least one association versus SNPs with no associ- 
ation. The result showed that the circadian set had a signifi- 
cantly higher ratio (OR=1.36, P=5e—29). Then, we 
calculated the total of unique traits associated with each 
of the variants, given that the SNP has at least one associ- 
ation. We used a Mann-Whitney U test to understand 
which set is represented by a higher level of traits per 
SNP. The circadian set was slightly more pleiotropic, and 
the result is supported by a significant P-value (P= 5.4e—3). 


Identifying Introgressed Circadian Variants with 
Evidence of Adaptive Introgression 


We sought to identify circadian variants that contain evi- 
dence of adaptive introgression. To achieve this, we col- 
lected adaptive introgression predictions from a method 
that applied various summary statistics on 1000 Genomes 
Project data (Racimo, Marnetto, and Huerta-Sanchez 
2017) and two sets of genomic regions that were measured 
for their likelihood to be under adaptive introgression 
by two machine learning methods: genomatnn and 
MaLAdapt. genomatnn is a convolutional neural network 
trained to identify adaptive introgression based on simula- 
tions (Gower et al. 2021). MalAdapt is a machine learning 
algorithm trained to find adaptive introgression based on 
simulations using an extra-trees classifier (Zhang et al. 
2023). Following the thresholds used in each paper, a re- 
gion is considered to be under adaptive introgression if 
the prediction value assigned to it meets a threshold of 
0.5 or 0.9, respectively. To find the variants of interest 
that fall into adaptive introgression regions, we intersected 
the set of introgressed circadian SNPs with the Racimo et al. 
(2015), genomatnn, and MalAdapt regions individually. 
The set of introgressed circadian variants contains variants 
inside circadian genes, in circadian promoter regions 
(5-kb upstream and 1-kb downstream of the TSS), and var- 
iants with regulatory function (cCREs) flanking circadian 
genes by 1 Mb. 


Supplementary Material 


Supplementary data are available at Genome Biology and 
Evolution online (http:/www.gbe.oxfordjournals.org/). 
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